Vibrational energy relaxation (VER) of isotopically labeled amide I modes in 
cytochrome c: Theoretical investigation of VER rates and pathways 

Hiroshi Fujisak^] 

Department of Chemistry, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215, USA and 
Institute of Physical and Theoretical Chemistry, J. W. Goethe University, 
Max-von-Laue-Str. 7, 60438 Frankfurt am Main, Germany 

John E. StrautQ 

Department of Chemistry, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215, USA and 
Department of Chemistry and Biochemistry, Montana State University, Bozeman, Montana 59111, USA 

(Dated: February 1, 2008) 

Using a time-dependent perturbation theory, vibrational energy relaxation (VER) of isotopically 
labeled amide I modes in cytochrome c solvated with water is investigated. Contributions to the 
VER are decomposed into two contributions from the protein and water. The VER pathways are 
visualized using radial and angular excitation functions for resonant normal modes. Key differences 
of VER among different amide I modes are demonstrated, leading to a detailed picture of the spatial 
anisotropy of the VER. The results support the experimental observation that amide I modes in 
proteins relax with sub picosecond timescales, while the relaxation mechanism turns out to be 
sensitive to the environment of the amide I mode. 



I. INTRODUCTION 



Amide I vibrational modes in proteins or peptides are localized around the peptide backbone and can be a sensitive 
probe of protein structure and dynamics. Mainly localized around the CO bond of the backbone with large oscillator 
strength, the amide I modes detected using infrared (IR) spectroscopy have been studied for a variety of proteins and 
peptides [H, [3, H, 0, Q ■ Recently 2D-1R spectroscopy [1, 0, H has been utilized to decipher (anharmonic) coupling 
between vibrational modes including the amide 1 mode [9|, llO| . The interpretation of the associated spectra is not 
necessarily simple, and will benefit from the development of increasingly accurate theoretical models [ill [l3| . 

Combining MD simulations with ah initio calculations, several theoretical groups have recently devised sophisticated 
methods to characterize the effect of inhomogeneity on the dephasing time, T2, of the amide I mode, using small 
peptides in water [3, [l^, [l^, [13, [3- Anharmonic frequency calculations have been performed for a peptide- like 
molecule, N-methylacetamide (NMA), in vacuum [l^, |20, [Hi, [13 and in water [l^. For the amide 1 mode of NMA 
in water, theoretical investigations of VER have been carried out using a quasi-classical method 26j and a quantum 
mechanical perturbation method [23|. However, while there have been many experimental studies [H, [24|, [2^ of 
vibrational energy relaxation (VER) of amide I modes in proteins, related to Ti, there are relatively few theoretical 
studies. 

In this work, we study the VER properties of amide I modes in a protein, cytochrome c, in water, and clarify the 
VER pathways using quantum- mechanical time-dependent perturbation theory (27j . This is a continuation of our 
previous work on VER of a CD stretching mode in cytochrome c [H, |2§| . The amide I modes studied were isotopically 
labeled according to IR experiment (s^, [3l|. We decompose the VER rate into two components: one from the protein, 
and the other from water (solvent). We identify the resonant modes that contribute most significantly to the VER 
rate, and introduce distribution-like functions for those modes to visualize the spatial anisotropy of VER as pathways 
in 3-dimensional space. 

This paper is organized as follows. In Sec. [Tll our time-dependent perturbation method is briefly described, and 
several methods to visualize the VER pathways in protein and water are proposed. In Sec lIIIi we examine the results 
of VER rates by comparing with the Maradudin-Fein formula, SASA, localization length, and experimental results. 
Furthermore, we mention the anisotropy of VER from excited amide I modes into protein as well as into water. In 
Sec. IIVI we summarize and discuss the future extension of our work. 
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II. METHODS 
A. VER rate formula 

We briefly summarize the quantum mechanical perturbation theory employed in this work. We assume that normal 
mode coordinates qa provide a good description of the system dynamics for the observables studied. We further 
assume that the potential energy function can be Taylor expanded up to the 3rd and 4th order anharmonic terms 
with respect to normal mode coordinates including only the relaxing mode qs [27j . 



H = Hs + HB-qsST + qlSg (1) 
2 



ns = 4 + Viqs) (2) 



^B-Ef + ^ (3) 
= ^Csa^iqaqp - {qaqp)) (4) 

Qp)) + ^CsSaqa (5) 

a, (3 a 

where Tis CHb) is the system (bath) Hamiltonian, and Csa/s (Cssas) are the 3rd (4th) order coupling terms. (This is 
related to the three- mode representation of a quartic force field |22l. |32||.) From the von Neuniann-Liouville equation, 
a reduced density matrix for the relaxing mode is derived using the time-dependent perturbation theory after tracing 
over the bath degrees of freedom. (A similar result has been derived from the path integral formulation of quantum 
mechanics by Okazaki and coworkers ^33,].) The Markov approximation is usually employed to derive a simplified 
Bloch-Redfield type equation after introducing the density of states for the bath. However, to describe the initial 
stage of the quantum dynamics, we do not invoke the Markov approximation. Using this approach, we are able to 
avoid assumptions related to the lifetimes of the vibrational modes of the bath — the so-called "linewidth problem" 
[28t used in the Markov approximation. 

When the relaxing mode is excited to the v = 1 state, the VER is described by the decay of the reduced density 
matrix element {ps)ii{t). We define the temporal VER rate as R{t) = d{ps)ii{t)/dt, which is approximately written 
as (23 



^(*) — 72 ["Cf^ut((l'g - UJg - LUfi) + C"^Mi(LJ5 + tJg + up) 

where ut{^) — sinr2t/f2, ujs is the anharnionicity-corrected system frequency, uja is the bath mode (harmonic) 
frequency, and the coefficients C"^^ , C"^ , C"'!, are derived from the nonlinear coupling constants Csap and Cssap 
[23] ■ R"^{t) is each component of the VER rate corresponding to modes a and /3. We note that the first term in 
Eq. ([6]) dominates in the formula because of the resonance condition, i.e. iit (f2) becomes large when (Dg — uja — up — 0. 

Because of non-Markovian properties, R{t) is not necessarily a constant. If the Markov assumption holds, R{t) 
becomes a rate constant. In the case of the amide I modes in cytochrome c, we have confirmed that a "VER rate" 
can be defined as a time average of R{t) after a certain transient time ttr 

i?^— i— f R{t)dt = Y,R-^ (7) 

using T ~ 1.0 ps and ttr — 0.5 ps. With this procedure, the averaged VER time is obtained as Ti ~ 1/R. The final 
rate formula is similar to the Maradudin-Fein (MF) formula [1^ [3J], and both formulas lead to a rate. An essential 
difference is that we avoid invoking the Markov approximation and introduciiig the phenomenological linewidth 
parameter. However, the MF formula has been widely employed in the literature [1^, HI], and it would be interesting 
to examine the usefulness of the MF formula based on our formula. Hence we introduce the MF rate as 

j^MF ^l^Y^ctt ^ ^ (8) 

^ Tp ~^ -^c- iop)^ 

where 6 is the linewidth parameter which will be discussed below. 
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B. Analysis of VER pathways 

As shown in Eq. (O, the VER rate R can be decomposed into components R"^ . This result aUows us to examine 
which modes (combination) contribute significantly to the VER rate. To develop a detailed picture of VER in a 
complex molecular system, we utilize the normal mode eigenvectors Uj^ which satisfy 



kl 



U^FkiUi' ^LulS^p (9) 



where a (k) represents the normal mode index (Cartesian coordinate index) and Fki is the mass-weighted Hessian 
matrix. 

We first introduce the localization length of the normal mode [s^, [IBl 



a |4 
k 



(10) 



which measures how a normal mode extends in a system in terms of the number of atoms. This measure (also called 
inverse participation ratio) is used to quantify the localization of the amide I modes. 

We are interested in evaluating contributions from protein and water to the value and mechanism of VER. We 
assign modes to water (protein) if the condition X^fcewator I'^fc 1^ > 0.5 (X^fceprotcin l^fc 1^ ^ 0.5) is satisfied. Here 
k £ water (protein) means that the summation is taken only for water (protein) degrees of freedom. Hence R"^ can 
be decomposed into the (1) protein-protein contribution (if a and (3 g protein), (2) protein- water contribution (if 
a{P) e protein and f3{a) e protein), and (3) water- water contribution (if a and /3 G water). 

We next introduce radial and angular excitation distribution functions for each normal mode a as 

5"(r) = J2 (11) 
fceAr 

h'^{0,4>) = \Ukf (12) 

where AT and Afi are the bins for the radial direction (whose median is r) and the bins for the cubic angle (whose 
median is 9, (/)), respectively (see Fig.[T]). The summation is only taken over these bins. We refer to <?"(?") and /i"(6', 0) 
as excitation distribution functions because they represent coarse-grained spatial information about the excitation of 
a normal mode. 

Note that these are different from the ordinary distribution functions (the absolute value does not have physical 
meaning), and are not directly related to experimental observables (such as neutron scattering functions). These 
functions are specifically introduced for visualizing the VER pathways. Other approaches are possible. For example, 
Dijkstra and Knoester visualized normal modes of /3-hairpin and /3-sheet peptides using a color map in their analysis 
of IR and 2D-IR spectra 

We define the average excitation distribution functions for resonant modes as 



ff(^) ^ ]^ E a'-ir) (13) 

a G res 
"'^ as res 

where a G res indicates that the summation is taken only for resonant modes, meaning those modes with R"^ > Rth- 
Nj-cs is the number of such resonant modes. We took Rth — 0.04 ps^^ for purposes of illustration. As indicated in 
Eq. when a system mode S is excited, the excess energy mainly fiows to resonant (bath) modes a and (3 via 
Fermi resonance {ujs — uja — ujp 0). We can explore the spatial content of such resonant modes using the excitation 
distribution functions. 

Sagnella and Straub introduced another way to visualize anisotropy of energy flow in a protein [36], and the 
differences are (1) we consider quantum mechanical energy flow whereas their calculation is classical, and (2) we 
employ the resonant normal modes whereas they used excess kinetic energies as the definition of the energy flow 
pathway. Okamoto and Nagaoka utilized an approach derived from hydrodynamics to visualize the anisotropic energy 
flow from a diatomic molecule in water [37| . though their work is also based on classical mechanics. Much similar to 
the present work is that by Mikami and Okazaki who devised a strategy for analyzing anisotropy of energy flow using 
quantum mechanics [38| . but the method is tailored for VER problems of a diatomic molecule. In this work, we take 
the naive approach described above to visualize resonant normal modes. 
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FIG. 1: The polar coordinates (r, 6, (j)) for an atom X used to calculate the radial and angular excitation functions of a normal 
mode. The origin of the coordinate system is the Cfs atom of the amide I mode. The sheet represents the amide plane. 

C. Numerical procedure 

We used CHARMM [s^ to construct and simulate the system consisting of cytochrome c and solvent water. 
CHARMM facilities were utilized to calculate dynamics, normal modes, solvent accessible surface area (SASA) etc., 
and the computational details have been presented elsewhere [13] ■ After equilibration, a 100 ps trajectory was 
generated and configurations were saved every 1 ps. From the 100 sample configurations, we calculated the VER rate 
using Eq. ([7]) and the results were averaged. 

Application of our formula starts from the decomposition of the system Hamiltonian using normal modes as in 
Eq. ID). However, computation of normal modes for the full system is prohibitive. Hence we take a reduced system 
strategy: We deleted all atoms except those in a spherical region r < Rc (see Fig.[T]) where Rc = 12 A [23] ■ For the 
reduced system, we computed the anharmonic coefficients Csap,Cssap using a finite difference approximation (28j 
and the normal mode frequencies tOa using instantaneous normal mode analysis (26l . [4]| (imaginary frequencies are 
neglected). The summation in Eq. ([7]) was taken only for the bath modes with uja > Wc, where we used Wc = 100 
cm~^ as the frequency cutoff. This is a reasonable approximation, as we have observed that the low frequency modes 
do not contribute significantly to the VER rate. 

In the CHARMM program, we isotopically labeled one CO bond in a residue as -'^•^C=^*0 [s^, [3l| and we examined 
four different cases as described below. It was expected that isotopical labeling should localize the amide I mode, but 
this is not necessarily the case as shown below. 

III. RESULTS AND DISCUSSIONS 
A. Protein and water contribution to VER rate 

We isotopically labeled four residues: 81st, 84th, 93rd, and 97th (see Fig. [2]). The first two belong to a loop region 
of cytochrome c whereas the latter two to a a helical region. By examing SASA (Table |l|, we see that the CO atoms 
of the 84th residue is most exposed to solvent whereas those of the 97th residue is most buried in the protein. 

In Table U we observe that the VER rates are ~ 2 ps~^ (Ti ~ 0.5 ps) for all the cases. This is similar to the VER 
rate of the amide I mode of N-methylacetamide in heavy water (~ 2.0 ps~^) [23| though the 97th residue has a slightly 
slower rate. It is interesting to note that the MF formula Eq. ([5]) gives rather "accurate" values with an appropriate 
linewidth parameter 5 = 3 cm^ ^ . This result seems to validate the use of the MF formula and to indicate that the 
MF formula is sufficient to describe VER. However, note that an appropriate value of the linewidth parameter 5 can 
not be determined in advance. 
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FIG. 2: Left: 81st and 84th residues in cytochrome c in a loop region. Right: 93rd and 97th residue in cytochrome c in a a 
helical region. The cartoon represents the protein, and the licorice the four residues. (The water molecules are excluded for 
simplicity.) VMD (Visual Molecular Dynamics) was used to generate these figures [53 |. 



TABLE I: Protein and water contributions to the VER rate for the amide I modes of four residues in cytochrome c in water. The 
total rate is R, and the protein-protein, protein- water, water-water contribution are denoted as Rpp,Rpnj and Rw-w, respectively. 
The value derived from the Maradudin-Fein formula is denoted as R^'^ . The units are in ps~^ for the VER rate. The inverse 
participation ratio (localization length) of the amide I modes and solvent accessible surface area (SASA) for the CO atoms are 
also shown. The units for SASA is A^. The number in the parenthesis represents the standard error of the mean. 



Residue 


R 


Rpp 


Rpw 






IPR 


SASA(CO) 


81st 


2.29 (0.03) 


0.95 (0.02) 


1.21 (0.02) 


0.13 (0.00) 


2.24 (0.03) 


9.54 (1.03) 


~ 80 


84th 


2.36 (0.03) 


0.99 (0.02) 


1.31 (0.03) 


0.06 (0.00) 


2.29 (0.04) 


8.26 (0.95) 


~ 100 


93rd 


2.27 (0.03) 


1.18 (0.02) 


1.05 (0.02) 


0.04 (0.00) 


2.24 (0.03) 


3.92 (0.03) 


~ 90 


97th 


1.75 (0.03) 


1.46 (0.03) 


0.30 (0.01) 


0.00 (0.00) 


1.82 (0.03) 


4.71 (0.13) 


~ 70 



Looking into the contributions from protein and water, we observe that the protein-water contribution is significantly 
less for the 97th residue compared to other residues (81st, 84th, and 93rd). This is because the water less surrounds 
the CO bond in this case (see Fig. [5]) as quantified by SASA for the CO bond (see Table There seems to be a 
correlation between SASA and RP"^ . However, as shown below, the water motion 8 A away from the CO bond can 
be involved in the VER processes. Hence the interpretation is not so simple. 

We also investigate this point in terms of the localization length Eq. (|10p for the amide I modes. The amide I 
mode is expected to be localized on a CO bond only but this is necessarily the case as shown in Table HI For the 
four residues we examined here, the amide I mode extends over from 4 to 10 atoms. For more delocalized modes, 
the protein-water or water-water contributions are expected to be large because there are more contact with water. 
However, this expectation is not validated in a strict sense as shown in Table |T1 We know that the Fermi resonance 
parameter |22l |29| can strictly describe the VER pathways in molecules, but for the present system it is difficult to 
interpret the result based on some chemical intuition. 

Note that the protein-protein contributions for the 81st, 84th, and 93rd residues are smaller than that for the 97th 
residue. However, due to the protein-water contribution, the VER rates of the former three amide I modes are faster 
than the latter. We conclude that VER of the former amide I modes is protein-mediated, whereas that of the latter 
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residue is water-mediated, which should have some experimental consequences as mentioned below. 

How do our calculated VER rates obtained in this theoretical study compare with experimentally derived values? 
Hochstrasser and coworkers obtained sub-picosecond VER time scales for amide I modes in several proteins [l^, 
[23 |. Zanni and coworkers also observed sub-picosecond VER time scales for isotopically labeled amide I modes in 
a membrane protein |3Q] . The experimentally derived VER times are comparable to our theoretical results (sub 
picosecond). However note that although the total VER rates are similar, the physical mechanism is sensitive to the 
specific structural environment of the amide I mode as discussed above (Table|T]). Using multidimensional spectroscopy 
0, d, [13, [13, [lOl, it should be possible to discriminate the protein and water contributions to VER separately in 
experiment. 

We have recently studied the solvent effects on heme cooling in myoglobin [4^ . We found that the mechanism for 
VER is dependent on resonance as well as structural anisotropy that modulates mode coupling. It will be interesting 
to investigate the solvent effects on VER of amide I modes by changing the nature of the solvent within the quantum 
mechanical perturbation theory (while the previous work 43] is based on classical mechanics). 

B. Visualization of VER pathways 

We have used the previously defined excitation functions Eqs. and to scrutinize the spatial anisotropy 
of the VER pathways. The radial excitation function in Fig. [3] shows each (protein or water) contribution to VER 
pathways along the radial direction. The peaks around r ~ 2 A correspond to the nearest neighbor intramolecular 
pathways involving the motion of the NH bond (see Fig.[T]). Since the density of the protein has a peak around r ~ 6 
A , the radial excitation function also has a peak or shoulder around at the same location. 

Looking at the water contributions, we found that the excitation functions have significant values at large distances 
(> 8 A), indicating that VER is mediated by distant water molecules. On the other hand, Mikami and Okazaki found 
that only the water molecules of the first solvation shell contribute to VER in the case of a solvated diatomic molecule 
[sst . We note that there are significant differences between the relaxation of an amide I mode in a protein and of a 
solvated diatomic molecule. For VER (or energy transfer) to occur, the overlap of resonant normal modes is needed 
[m . In the case of a diatomic molecule, the stretching mode is, by definition, localized around the bond, and mainly 
couples to the neighboring water molecules. However, in the case of amide I modes in solvated proteins, the amide I 
modes are delocalized as shown in Table [H Many resonant modes couple to such amide I modes, which can contain 
delocalized (collective) water motions. Further theoretical and experimental studies are needed to clarify the role of 
such delocalized modes on VER in proteins. 

Finally we examine the angular excitation functions in Figs. [4] and [5] to understand the VER pathways along the 
angular directions. We observe that the VER of the amide I modes is spatially anisotropic. For the protein (Fig.[3]), 
the resonant modes are mainly localized around (/) = or ±7r, which indicates that VER occurs in the peptide 
backbone plane (see Fig. [J). For the 81st and 84th residues, belonging to a loop structure, the value around (/) ~ 
enhances because the amide plane is planer and the intramolecular contribution comes from such a direction. For the 
93rd and 97th residues, belonging to an a helical structure, the value around (j) diminishes because the amide I 
plane is bended. 

On the other hand, for the water (Fig. [5]), the angular excitation function is rather broad, and these distributions 
are similar to the conventional distribution function of water around the CO bond (not shown here). However, 
the excitation functions are less uniform compared to the conventional distribution function of water, implying the 
anisotropy of VER into water. While these excitation functions can not be directly measured, it is of significant 
interest to examine the anisotropy of energy flow in proteins using experimental methods such as the multidimensional 
spectroscopy 0, S, 0, [M [131 ■ 

Recently Dlott and coworkers [31 and Hamm and coworkers [i^l devised a new experimental technique to clarify 
energy transfer pathways in molecular systems. This method may be also combined with our analysis to clarify the 
water contribution of VER in a protein. 

IV. CONCLUDING REMARKS 

Using the time-dependent perturbation formula we developed in our previous paper, we have investigated the 
vibrational energy relaxation (VER) of amide I modes in cytochrome c solvated with water. We observed that the 
VER rates are sub-picosecond for the four residues we examined here, which are in accord with previous experiment for 
other proteins. We have decomposed VER into separate contributions from protein and water, and further projected 
the resonant modes onto radial and angular excitation functions. Although the total VER rate is similar, the detailed 
mechanism of VER is different; there are protein- mediated pathways and water-mediated pathways depending on 
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FIG. 3; Radial excitation functions for the resonant normal modes of protein and water for the (a) 81st, (b) 84th, (c) 93rd, 
and (d) 97th residues, represented in arbitrary units. 



which residue is excited. There is anisotropy of VER for the water-mediated pathways as well as for the protein- 
mediated pathways, implying experimental consequences by multidimensional spectroscopy or other means. 

Anharmonicity is a key issue in VER calculations. We note that our simulated dynamics are based on the CHARMM 
force field, in which vibrational anharmonicity has been underestimated. There is some evidence that the anharmonic 
coupling calculated using the force field can be comparable to that derived from ab initio calculations (Fujisaki, 
H.; Yagi, K.; Hirao, K.; Straub, J.E. unpublished). However, without such a direct comparison, it is unclear how 
accurate the results of our force field calculations are. Nevertheless, the use of empirical energy functions is currently 
the only feasible way to characterize VER in large molecules. To expand on these studies, we must develop a new 
methodology to investigate the quantum dynamics of the amide I mode with ab initio potential surfaces. As first steps, 
we investigated N-methylacetamide in vacuum using the VCI method and in water cluster using the perturbation 
method (Zhang, Y.; Fujisaki, H.; Straub, J.E. unpublished) on ab initio potentials. How to extend these studies to 
larger systems using multiresolution methods [i^lJ^I or QM/MM methods [i^, H^ll^l would be our focus in the near 
future. 

Another important and related topic which should be pursued by such ab initio methods is the effect of polarization 
[13, [iM IitI- Tsj . Morita and Kato clarified that polarization is important for VER of a stretching mode in azide 



ion in water [51[, and this might be also the case for VER of amide I modes in a protein. In particular, it is important 
to investigate some portion of a protein where the internal electric field is large and the polarization effect is expected 
to be enhanced. 
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FIG. 4: Angular excitation functions for the resonant normal modes of the protein for the (a) 81st, (b) 84th, (c) 93rd, and (d) 
97th residues, represented in arbitrary units. 
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